Kramers Theory for Conformational Transitions of Macromolecules 
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We consider the application of Kramers theory to the microscopic calculation of rates of confor- 
mational transitions of macromolecules. The main difficulty in such an approach is to locate the 
transition state in a huge configuration space. We present a method which identifies the transition 
state along the most probable reaction pathway. It is then possible to microscopically compute the 
activation energy, the damping coefficient, the eigenfrequencies at the transition state and obtain 
the rate, without any a priori choice of a reaction coordinate. Our theoretical results are tested 
against the results of Molecular Dynamics simulations for transitions in a 2-dimensional double well 
and for the cis-trans isomerization of a linear molecule. 

PACS numbers: 



The kinetics of conformational changes of macro- 
molecules is believed to provide important information 
about the underlying mechanisms involved in such re- 
actions. In such a context, rates are the fundamental 
observables. Not only they provide direct tests for the- 
oretical calculations, but they also encode information 
about the structure of the important reaction pathways. 
For example, by (^-value analysis it is possible to iden- 
tify the residues which are structured at the transition 
state 

Kramers theory and its multidimensional generaliza- 
tion offer a scheme to compute the transition rates for 
bistable molecular systems. In such a formalism 0, 
the transition rate of a particle in an external poten- 
tial U {x) in N dimensions, from the meta-stable state h 
across the saddle-point o can be written, in the strong 
friction regime, as 
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where 7 is the friction coefhcient, lJs is the angular fre- 
quency of the single unstable mode at the saddle-point 
and and w° are the stable frequencies in h and in o, 
respectively. The ratio ujsll is often called the (adimen- 
sional) damping factor. It is responsible for lowering the 
actual rate from the theoretical upper limit, , given 
by transition state theory. 

Kramers theory has been successfully applied to more 
complicated chemical reactions involving macromolecules 
[sl, such as two-state protein folding. In this context, 
Eq-© is used as a phenomenological prescription in 
which X is a set of reaction coordinates, U{x) is the cor- 
responding potential of mean force (free energy) and 7 is 
the effective friction at the transition state. 

In the present work we adopt a different strategy: The 
configuration x specifies the microscopic degrees of free- 



dom of the molecule (e.g. the atom or residue coordi- 
nates) and U{x) is the interaction potential with implicit 
solvent [sj. The main difficulty in such an approach to 
chemical reaction rates, is that it requires to know the 
location of the saddle-point state a; = o in a large config- 
uration space. This information is needed to determine 
the activation energy, the damping factor and the eigen- 
frequencies, which are in turn needed to estimate the rate 
constant. In practice, the identification of the transition 
state in molecular reactions represents a very challenging 
task and Kramers formula cannot be directly applied. 

The key point of this work is to show that the problem 
of finding the transition state in two-state conformational 
reactions of macro-molecules can be efficiently solved us- 
ing the recently developed Dominant Reaction Pathways 
formalism 0, 0, H, @] . This formulation of the stochastic 
dynamics leads to an impressive computational simpli- 
fication of the problem of finding the most irnportant 
reaction pathways in high dimensional systems [6|. The 
reason is that the dominant reaction pathway is sampled 
at equally-spaced displacement steps, rather than using 
constant time steps. In thermally activated reactions, 
due to the decoupling of time scales, the difference be- 
tween these two samplings is huge. In particular, for the 
folding of a polypeptide chain, the number of displace- 
ment discretization steps is of order 30 — 100 0, 0] . This 
number should be compared with the order 10^^ steps 
which would be required to describe the same reaction 
using constant time steps. In our recent work 0,0, we 
have shown that it is possible to determine the most sta- 
tistically important reaction pathways in conformational 
transitions of amino-acid chains. In this Letter, we show 
how to perform the atomistic calculation of Kramers re- 
action rates in macromolecular transitions by using the 
most probable paths, with available computers. 

We begin by briefly reviewing the Dominant Reac- 
tion Pathways approach (for a detailed and self-contained 
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introduction, see [9|]). For the sake of simpUcity, we 
present all formulas in the case of the diffusion in a one- 
dimensional external potential. The generalization to the 
multidimensional case is straightforward. Let us there- 
fore consider the overdamped Langevin dynamics of a 
system with coordinate x in a potential U{x): 
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where D ~ ksT/"/ is the diffusion coefhcient, and ksT 
is the thermal energy. r]{t) is a Gaussian noise with zero 
average and correlation given by {ri{t)ri{t' )) — 2Dd{t—t'). 
Note that in the original Langevin equation an inertial 
term, mx, appears. However, in the case of the dynam- 
ics of biopolymers in water, the effect of such a term 
can be neglected for time scales larger than a fraction of 
picosecond [l^. As is well known, the stochastic differ- 
ential equation (U) generates a probability distribution 
P{x, t) which obeys the Fokker-Plank Equation (FPE) 
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In the study of noise-driven reactions, we are interested 
in transitions between two meta-stable states a and b. 
The starting point of the Dominant Reaction Pathways 
approach is to express the solution of the FPE, subject to 
the boundary conditions a;(0) = b and x{t) = a in terms 
of a path integral: 



P(6,0;a,t) 
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where 5'ef[2:] — J dr (^x^ /AD + Vej f[x]) is an effec- 
tive action and the effective potential reads Veff{x) = 
D/{2kBTf [{VU{x)f - 2kBTV^U{x)] . 

The most probable paths in configuration space con- 
tributing to (ID) are called the dominant reaction path- 
ways (DRP). They are those for which the exponential 
weight e~^^ff is maximum, hence for which the effective 
action Seff is minimum. In our recent works 0,0,0] we 
have shown that the DRP can be rigorously obtained by 
minimizing the effective Hamilton- Jacobi (HJ) functional 



SHA[x];b,a)^ di^D-^Veff[x{l)]~Veff{a)), (5) 



where dl = y^jjdx^ is a measure of the elementary 
displacement along the reaction path. We note that the 
HJ functional does not depend on the diffusion coefficient 
D. Hence, the DRPs are not sensitive to the choice of 
the friction coefficient. This could be seen already at the 
level of the FPE ([3|) , in which the choice of the diffusion 
constant just sets the time scale. 

The transition state can be identified as the configu- 
ration for which the potential energy is maximum along 




FIG. 1: Contour plot of the two-dimensional energy landscape 
used to test the method. 



the DRP, and is a saddle point on the full potential en- 
ergy landscape. Once the transition state Xts has been 
found, one can easily compute the eigenfrequencies and 
damping factor required to compute the rate from eq. 

With this knowledge, it is possible to estimate the 
friction coefficient in the transition state, which enters 
in ll]), by means of short MD simulations with explicit 
solvent. 

Before showing examples of applications of the present 
method to specific systems, we first show that the DRP 
contains much more information about the kinetics of the 
reaction than encoded in the Kramers rate formula. In 
fact, it is possible to use this method to obtain micro- 
scopic parameter-free predictions not only for the rate, 
but also for the probability of performing an arbitrary 
number q of transitions between b and a, in a given 
time_interval t. We start from the work of Caroli et 
al. 11], who studied the diffusion in an asymmetric one- 
dimensional double-well, by estimating the path integral 
dH) in the so-called dilute instanton gas approximation. 
They obtained the expression for the transition probabil- 
ity form 6 to a in a time t, with unconstrained number 
of barrier crossings: 



P(6,0;a,t) = 
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ki = Dy/\U"{o)\ U"{i){2TTkBT)-^e~"^'''B"^''' . (7) 

Note that the relaxation to the equilibrium distribution 
is controlled by the rate {ka + kt), which is precisely 
Kramers rate. A generalization of this result to a fixed 
number q of barrier crossing in time t is made possible by 
exploiting the peculiar form of the path-integral solution, 
which allows to find the analytical solution for a generic 
double well. Here we show only the result for the sym- 
metric double-well, which takes the particularly simple 
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LJs/^Tf,q, where and r^q are 



FIG. 2; Comparison between Dominant Reaction Pathway 
prediction (dashed line) and MD simulations (solid line) for 
the P''')(f), with g = 1, 3, 5, 7, in the two-dimensional double 
well of Fig. [U 



form P(9)(t) oc {kbtf e-'"'* /q\. It is easy to check that 
after summing over all possible numbers of crossing, one 
recovers the correct unconditioned probability ^ with 
Kramers time = (2/cti)~^ now appearing at the expo- 
nent. These results can be generalized to bistable sys- 
tems in higher dimensions, by replacing the expression 
in Eq.([7]) for fcf, with 
defined as in Eq.([T]). 

We now illustrate how the present method is imple- 
mented in practice and we compare it with Molecular 
Dynamics (MD) simulations results, focusing on the gen- 
eralized transition probabilities P'-'^' (t) , since they pro- 
vide a very stringent test for the approach we have de- 
veloped. 

Consider the simple case of the 2-dimensional bistable 
potential U{x,y)/kBT = A{x^ - if + By^, plotted in 
Fig. [T] with A = 4 and B = 15. In Fig. [2] we report 
the sampled histograms of P'^)(t) up to g = 7, obtained 
from several MD simulations, and we compare it with the 
corresponding theoretical curves obtained by the present 
analysis. In this case the dominant trajectory is unique 
and coincides with the straight line connecting the min- 
ima of U{x, y). 

On the other hand, the construction of the histograms 
requires a prescription to record crossing events. We de- 
cided to set two thresholds, one on the left and one on 
the right of the barrier top, at x = ± 0.4, respectively. 
A crossing event is marked only when the two thresholds 
are crossed in sequence. The histogram of the transi- 
tions is updated each time the system is closer to a than 
a given distance S. The domain width S and the inte- 
gration time-step have been set respectively to 0.02 and 
3 X lO-'^. 

The agreement between the analytic calculation based 
on the DRP and the histograms obtained from the MD 
simulation is excellent and holds for all the conditional 



probabilities, P^^-*(t), P^'^\t). We emphasize again the 
fact that the theoretical curves are parameter-free and 
therefore the agreement represents a compelling evidence, 
supporting the validity of the method. 

Let us now consider the application of the Dominant 
Reaction Pathways method to the kinetics of conforma- 
tional transitions of a molecular system. It is important 
to stress the fact that for sufficiently complex molecules, 
the position of the transition state in a conformational 
transition cannot in general be inferred a priori. This 
limitation prevents the direct application of Kramers the- 
ory. The aim of the present analysis is to check if the 
present approach is able to find the correct saddle-point, 
and predict the P^'^^{t) (and therefore Kramers rate). To 
this end, we consider the cis-trans isomerization of a toy 
molecule, in which the exact location of the transition 
state is known by construction. In particular, we used 
a model for a linear molecule composed of 8 interaction 
centers, with fixed bond lengths of 1 A, masses of 10 
a.m.u., and torsional potentials acting on the 5 dihedral 
angles (fid, of the form U — Cd[l + cos{nd(l)d — 0^)] • The 
interaction parameters are C3 = 3kJ/mol, = 2 and 
03 = 90° for the central dihedral (j>3 and Cd = lOkJ/mol, 
rid = I and (j>'^ — for the remaining four dihedral angles. 
The plane angles between every three consecutive atoms 
are kept fixed at 90°. This potential has two minima, 
located at the cis and trans conformation relative to the 
central dihedral, respectively. 

We computed the DRP by minimizing numerically the 
HJ action by means of a simulated annealing algorithm, 
starting from several randomly generated arbitrary paths 
connecting the two states. We then looked for the max- 
imum of the potential energy along the reaction path to 
identify xts ■ The correct location of the saddle-point took 
only few minutes of CPU time and gave the same result, 
regardless of the random starting path used. 

The next step consists in computing the equilibrium 
constant 1 /rgq which enters the expression of the P(«)(i) 
and of the Kramers rate formula. If the only active de- 
grees of freedom are the torsional angles (bond length 
and angles are kept fixed), the Jacobian of the coordi- 
nate transformation from Cartesian to internal ones does 
not depend on the internal coordinates themselves [l^ . 
This fact can be used to simplify the expression for the 
equilibrium rate constant for a generic molecule Q to the 
following: 



1 _ 1 \M-^/^WJ{o)\ 
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(8) 



Here TV is the number of the internal degrees of freedom, 
M is the diagonal matrix of the atom masses, and x and (f) 
indicate the Cartesian and internal coordinates (dihedral 
angles in the present case), respectively. The (positive) 
eigenvalues of the Hessian matrix Hij = ^aFi'^. evaluated 
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FIG. 3: Comparison between the P''{t) functions for the cis- 
trans isomerization of the toy molecule, obtained from MD 
simulations (solid lines) and from the DRP (dashed lines). 



at the saddle point and in the starting weU are denoted 
A° and A**, respectively. The damping factor F can be 
estimated H as F = --A= ^^^IZ^^^i-^ , where A is 

the negative eigenvalue of the D(o)H(o) matrix, D(o) 
is the diffusion tensor in internal coordinates and, A° 
is the negative eigenvalue of the Hessian matrix at the 
transition state. The Kramers rate has to be estimated 
in the moderate-strong friction regime, where the adi- 
mensional damping factor F ~ 0.5. The computation of 

-1 



(2F)~^ + 1 - (2F)" 



T„,} , resulted in a rate of 



0.034 ps-i. 

As before, we are interested in assessing the reliability 
of this approach, by comparing the prediction of the gen- 
eralized transition probabilities with the results of MD 
simulations. We integrated the full Langevin equation 
using a diffusion coefficient of O.IA^ ps^^ and integra- 
tion time step of 0.001 ps, at a temperature of 300 K. 
The histograms obtained from 10^° integration time-step 
are shown in Fig. [31 along with the prediction obtained 
from the present method. The thresholds and domain 
width S for the cosine of the torsional angles have been 
chosen to be ±0.5 and 0.05, respectively. The agreement 
between the theoretical predictions and the results of nu- 
merical simulations is quite good, given the approxima- 
tions involved in the estimate of the rate constant. The 
calculation of the P^'i\t) starting from the DRP took 
only minutes of CPU time. On the other hand, approx- 



imatively 150 CPU hours were required to reconstruct 
the same curves from histograms obtained from MD sim- 
ulations. Note that such a gain was observed also for 
reactions in more realistic systems, using empiric atom- 
istic force fields Q- 

In conclusion, in this work we have developed a 
parameter-free method to compute Kramers rate at the 
microscopic level, in high-dimensional systems exhibit- 
ing two-state kinetics. Our theoretical expressions for 
the generalized transition probabilities p(')(<) have been 
found in excellent agreement with the results of MD sim- 
ulations performed in a two-dimensional double well and 
in the cis-trans isomerization of a simple molecule. This 
approach is very accurate and leads to a huge computa- 
tional gain with respect to MD simulations, and makes 
it possible to perform calculations of rates for systems in 
which MD techniques are not presently feasible. 

Computations were performed on the HPC facility at 
the Department of Physics of the Trento University and 
at the Frankfurt Center for Scientific Computing, whose 
support we gratefully acknowledge. We thank W. Eaton 
and A.Szabo for important discussions. 
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